Thermodynamic stability of small-world oscillator networks: 
A case study of proteins 
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We study vibrational thermodynamic stability of small- world oscillator networks, by relating the 
average mean-square displacement S of oscillators to the eigenvalue spectrum of the Laplacian 
matrix of networks. We show that the cross-links suppress S effectively and there exist two phases 
on the small- world networks: 1) an unstable phase: when p <C 1/N, S ~ N; 2) a stable phase: when 
p S> 1/N, S ~ p _1 , i.e., S/N ~ E~r . Here, p is the parameter of small-world, N is the number of 
oscillators, and E cr = pN is the number of cross-links. The results are exemplified by various real 
protein structures that follow the same scaling behavior S/N ~ E~ r of the stable phase. We also 
show that it is the "small-world" property that plays the key role in the thermodynamic stability 
and is responsible for the universal scaling S/N ~ E^r, regardless of the model details. 



PACS numbers: 87.14.E-, 05.40.-a, 89.75.-k 
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Vibrational dynamics has been widely used to study 
thermodynamic properties of various structures in solid 
state physics and/or other disciplines Since the 

structure is considered as a primary factor responsible 
for physical properties, to keep the underlying structure 
thermally stable is of primary important for systems to 
function properly. For example, proteins, comprising 
an extremely heterogeneous class of biological macro- 
molecules, must be stable enough against thermal fluc- 
tuations and/or external perturbations so as to maintain 
their native structures and to function correctly @, 
Therefore, a natural and important question is often 
asked: what is the structure effect on thermodynamic 
stability? In this paper, we study the stability of small- 
world structures [J|, which is then exemplified by pro- 
teins. 

The dynamics of N coupled oscillators on the network 
in contact with the external heat reservoir can be ex- 
pressed as: 



Mq = -aLq - Tq + £ 



(1) 



where q = [qi, q2, 9at] t , denotes the oscillator's dis- 
placements from the equilibrium positions. My = rriiSij 
is the mass matrix, where denotes the mass of the zth 
oscillator, a is the spring constant. Ly = £y Yl m ~ 
Aij is the Laplacian matrix and Ay is the adjacency ma- 
trix of the network, where Ay = 1 if i and j are con- 
nected and = otherwise. Ty = Yj<5y is the dis- 
sipation matrix where ji is the dissipation coefficient of 
the ith oscillator influenced by the heat reservoir. Vector 
£ = £2, £n] T denotes the thermal fluctuation with 
zero mean and variance (£i (£)£/(*')) = 2fcsTTy5(f — t'), 
which is the usual fluctuation-dissipation relation. 
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The harmonic potential we adopt looks very simple 
but can capture the main features of the system. For ex- 
ample, Tirion Q demonstrates that a single-parameter 
harmonic potential can reproduce vibrational properties 
of the real macromolecular system very well. Thereafter, 
the Gaussian network model (GNM) [6[ has been widely 
used in protein research and yields results in good agree- 
ment with experiments. In the GNM model, the interac- 
tions are considered as homogeneous harmonic springs, 
which is in analogy with the elasticity theory of random 
polymer networks 0, § • 

The correlation matrix of oscillator displacements at 
the steady state for Eq. (1) can be easily obtained (see 
Appendix |X| : 
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where matrix G{±iuS) = (±ioj) 2 AI + (±iw)F + <jL. Since 
G{iuj) — G(—ioj) — 2iojT and G(0) = <jL, one can elimi- 
nate T in the above integral and obtain: 
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id 



(3) 



where L' denotes the pseudo-inverse of L. It excludes 
zero mode which correspondes to the translational in- 
variance of the system, and is the inverse of L in the 
subspace orthogonal to the zero mode: 
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where A a are the non-zero eigenvalues, and ip a j denote 
the corresponding normalized eigenvectors of L. There- 
fore, we can obtain the average mean-square displace- 
ment straightforwardly: 
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FIG. 1: (al) The numerical result of the correlation of pairwise oscillator displacements on a ring chain N = 64. (a2) The 
numerical result of the correlation for adding link (11, 24) on (al). (a3) is the numerical calculation of correlation for adding link 
(21, 40) on (a2). (a4) is the correlation for adding link (33, 64) on (a3). (b2), (b3), (b4) are the perturbed calculation results 
using Eq. (8), corresponding to their counterparts in the upper panel, (bl) shows the decreasing of 5* with adding cross-links. 
The red ones are calculated directly from numerical diagonalization and the black ones are the perturbed calculation using Eq. 
(8). All the results validate our analytic results. 



This formula relates the dynamic vibration property S to 
the static structure property — the eigenvalue spectrum 
of Laplacian matrix L. When the average mean-square 
displacement S reaches the square of the typical spac- 
ing between oscillators, the structure encounters large 
vibrations and becomes unstable. Thus, small value of 
S means stable, while large value means unstable. It 
is clear that S has a trivial dependance on T and a so 
that lower temperature or larger spring constant indi- 
cates more thermal stability. Therefore, to study the 
structure effect on S, we take ksT/a = 1 in the follow- 
ing, without loss of generality. 

Based on Eq. ([3]) and ([5]), we can use exact numerical 
diagonalizations of the Laplacian matrix L to study the 
structure effect on thermodynamic stability. 

In fact, the stability can be also studied by pertur- 
bation analysis, through which we find that the cross- 
links can suppress the thermodynamic instability, i.e., 
decrease S effectively. After structure changes, the new 
Laplacian matrix is constructed as V = L — A, where 
A denotes the perturbation. The new correlation matrix 
C can be written as 



C = (I- CA) -1 C = C + GAG + CACAC ■ 



(0) 



This is a standard algebraic algebra treatment of matrix 
perturbation, which can be simply regarded as Taylor 
series. In fact, it is similar to the Dyson equation in 



Feynman-Dyson perturbation theory Q . For the case of 
adding a link between nodes i and j, the perturbation 
matrix A can be expressed as 

A mn — dmidnj ~\~ fimj^ni $mi$ni ^mj^nj ■ (7) 

Substituting the expression of A into C' , one obtains 



l + Ri 



(8) 



where Rij = Ci, 



C 



J 3 



Therefore, the new average mean-square displacement is 



S' 



— tr C = S - ^-*k=i(Cik - Cjk) 



S- 



1 



JV 



(9) 



The second term in Eq. ([9]) is positive so that the value 
of new S' is always smaller than S. In other words, the 
cross-links always decrease S so as to increase the ther- 
modynamic stability of the system. A specific case is 
studied and illustrated in Fig. [T] 

In the following study, we choose a typical model to 
construct the small-world structure 10 . We first con- 
sider TV oscillators (which might be an atom, a molecule, 
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FIG. 2: (color online), (a) S versus network size N. The 
straight line for p = indicates the diverging behavior, S ~ 
JV, of the ID ring chain. Even small nonzero value of p can 
suppress the diverging behavior of S and make it saturated 
to a finite value in the large size limit, (b) For each N, S 
decreases as p increases. A scaling behavior S ~ 1/p emerges 
in the large size limit. All the results indicate that the cross- 
links boost the thermodynamic stability effectively. 



or other module structure, depending on the system stud- 
ied) on a one-dimensional(lD) ring chain, i.e., with peri- 
odic boundary conditions. Each oscillator is connected to 
its nearest-neighbors. Then, we add a cross-link to each 
oscillator with probability p, which connects to another 
non-neighboring oscillator randomly Thus, E cr = pN is 
the number of cross-links. When p = 0, the structure 
reduces to the ID ring chain. In all cases studied below, 
each data point is obtained by averaging over 50 different 
network configurations for a given p and N. 

Figure HJa) illustrates S versus the system size N for 
different values of p in double logarithmic scale. The 
p = case, corresponding to the ID ring structure, shows 
the power-law divergence, S ~ N. It indicates that no 
thermodynamically stable solid exists at finite tempera- 
ture in ID. Indeed, when the average mean-square dis- 
placement S exceeds the square of the typical spacing 
between oscillators, the structure behaves like a liquid 
rather than a solid, and the crystalline order makes no 
sense anymore. This behavior is also reported in [ll|, E2| • 
For the case of p ^ 0, even of small value, as N increases, 
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FIG. 3: (color online). Scaling plot of the average mean- 
square displacement S in small- world networks for various size 
N and probability p. All data are from Fig. [2] 13] and they 
collapse into one single line very well. It shows two phases 
clearly: one is the regime with slope —1, where p 3> 1/iV, 
indicating the non-divergent stable behavior S/N ~ E^ r = 
(pN)" 1 , i.e., S ~ 1/p. Another one is the horizontal regime, 
indicating the diverging unstable behavior, S ~ N. The red 
vertical dashed line is used for eye-guiding to separate the two 
phases. 



S is saturated to a finite value rapidly. Moreover, from 
Fig. [2fb) , we can see that the larger p, the smaller S 
and in the large N limit, a scaling S ~ p _1 emerges. All 
the above results indicate that the cross-links suppress 
the average mean-square displacement S effectively and 
make it convergent in the thermodynamic limit. 

To eliminate the finite size effect, all the data from 
Fig. [2] [l3| are re-scaled and the results are illustrated 
in Fig. [3] It is found that all data points collapse into 
one single line very well and two distinct phases emerge. 
For p <C X/N, there is a horizontal regime, S/N ~ 
const. It indicates an unstable phase: when the num- 
ber of cross-links is smaller than one, S diverges with 
N. For p ^> 1/N, there is a regime with slope — 1, 
S/N - E-} = {pN)- 1 , which indicates a convergent 
stable phase. In other words, when the number of cross- 
links is much larger than one, S approaches a finite value 
at large N and scales as 

Since the average mean-square displacement S is re- 
lated to the spectral properties of the Laplacian matrix 
L, we can understand the scaling behavior of S in terms 
of its eigenvalue spectrum p(X). For large size N, Eq. ([5|) 
can be expressed as fl2j |: 



S = 



A 



dX, 



(10) 



from which we can easily see that the density of small 
A dominates the behavior of S. For the case with- 
out cross-links, the system reduces to a ID ring chain, 



4 



where p(X) ~ A -1 / 2 and A ~ N~ 2 for small A. Thus, 
S ~ J A~ 3//2 e?A ~ iV. For the case with cross-links, 
following the heuristic argument in [3], we can con- 
sider that the ring chain is divided into several quasi- 
linear segments of length /, and the probability of length 
I is exponentially small, e~ pl . Each segment I con- 
tributes to small eigenvalues of the order of Z~ 2 . Sum- 
ming over lengths with the exponential weight, we ob- 
tain S = J ~ l -^dl = 1(1 - e"^). When 

pN <C 1, S ~ N; while pN > 1, S ~ 1/p, which is 
exactly what our numerical results show in Fig. [2] and 
[3j Although the argument above is not rigorous and ap- 
plies only when p is smaller than one, it gives us quite 
good understanding of the scaling behavior of S. When 
p is larger than one, the model we used is more like an 
Erdos-Renyi model, which is also to be demonstrated to 
follow the same scaling of the stable regime at the end of 
this paper. 

The small-world structure we used above is well stud- 
ied [l(|. Usin g re normalization group method, the au- 
thors in Ref. [lOj showed that this model undergoes a 
transition between regular lattice and random one at in- 
termediate characteristic size N c ~ p~ x . In other words, 
the phase transition has a critical point p c = in the 
thermodynamical limit when N — > oo. For finite size N , 
the diameter I scales linearly with N for jV c N as it 
is in ID ring chain, while / ~ In N for N 3> N c , where it 
exhibits "small-world" property. Our results about un- 
stable and stable phases are consistent with their findings 
that the unstable regime corresponds to the ID case and 
the stable phase corresponds to the "small-world" case. 

As an illustrative example, the thermodynamic stabil- 
ity is further tested on real protein data. We revisit the 
proteins used in Ref. Q, which differ in functions and 
structures, with a wide size scale ranging from 100 to 
3600 residues. All the structure data are downloaded 
from the Protein Data Bank (PDB) [15| and the number 
of all residue pairs is counted within a customary cut- 
off 7.0 A. After eliminating the connectivity number of 
the primary structure of protein from the counted num- 
ber, we obtain the number of cross- links, E cr , for each 
protein. The mean square displacement of C Q atoms is 
characterized by S-factor [l6|, also called Debye- Waller 
or temperature factor: Bi = 8n 2 (q 2 )/3, where i is the 
index of amino acid residue. It is experimentally mea- 
sured via x-ray crystallography, and also can be down- 
load from the PDB. The average B- factor is calculated 
over all C Q atoms for each protein, B = Yli=i Bi/N. No- 
tice that at above theoretical analysis, ksT/a is set to 
be 1 for convenience, which is not always true. The value 
of ksT/a varies among different proteins. Thus, we use 
the estimated data of k B T/a Q to obtain the normal- 
ized average B-factor, B' = B/(k B T/a). Note that B' 
is analogous to S, defined in Eq. (fit)]) . All the details of 
these proteins are listed in Table Q] in Appendix [B] 

The thermodynamic stability is crucial to keep the na- 
tive structure of protein for right function. Moreover the 
structure of protein is also found to have "small-world" 
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FIG. 4: (color online), (a) Log-log plot of normalized average 
B-factor with various number of cross-links for real protein 
data. It exhibits clearly a power-law behavior, B' /N ~ E~ r a . 
The dashed line indicates the best-fit of the power-law, with 
exponent a = —0.92 ± 0.01. (b) The average mean-square 
displacement S versus the parameter E cr /N = p r (N — 1)/2— 1 
in Erdos-Renyi model. The dashed line indicates the scaling. 



property (l7| . i.e., I ~ IniV. It is intuitive for us to ex- 
pect that nature selection forces proteins evolving into 
the stable phase in Fig. [31 which implies B' /N ~ E^ 1 . 
Figure 0]^a) verifies our expectation drawn from the ar- 
gument of stability analysis. In fact, we obtain a clear 
power-law scaling: 



B'/N ~ E- c 



0.92 ±0.01, 



(11) 



which is quite close to 1. This scaling reveals the univer- 
sal behavior shared by various different proteins, regard- 
less of their sources or functions. It implies an underly- 
ing general mechanism that nature selects proteins with 
thermodynamic stability constraints. 

Although protein has more complex structures with 
high modularity (domains), "small- world" property cap- 
tures its main feature. Thus, "small-world" might play 
a key role in the thermodynamic stability of structures 
and be responsible for the scaling in the stable regime. 
To validate our conjecture, we further study the thermo- 
dynamic stability in Erdos-Renyi (ER) random network 
model in the following. 

ER model [l8( has N nodes and every pair of nodes 
is connected with probability p r . The average degree 
(k) = p r (N — 1). There are several phases in this 
model depending on different threshold p r : when (k) = 
p r (N — 1) > 3.5 [l8j], the diameter of the graph equals 
the diameter of the giant cluster, and is proportional to 
IniV, i.e., "small- world" property. Thus, it is straight- 
forward to expect that ER model might share the same 
behavior S/N ~ E^ 1 . The numerical result is illustrated 
in Fig. 4(b). As we point out above, ER model has 
"small- world" property pi} when p r (N — 1) > 3.5. Cor- 
respondingly, whenp r (JV-l)/2-l > 0.75, S ~ N/E cr = 



5 



io- 1 



10"' 



10' 



CO 



10' 



10"' 



■ small-world model 

• ER random model 

* real proteins 



\ 



10 10 10 10 10 10 10 10 



FIG. 5: (color online). S/N versus E cr for three different 
networks. It shows that the "small-world" property is re- 
sponsible for the universal scaling S/N ~ E^ r in the stable 
regime, regardless of the model details. Note that for proteins, 
S — B'/(8iv 2 ), where the factor 3 is removed since B-factor 
is measured in 3-dimension. 



N/( Pr N(N - l)/2 — N) = [ Pr (N - l)/2 - l]" 1 as shown 
in Fig. 4(b). 

For convenience of comparison, we plot the data of 
three cases together in Fig. 5. It clearly shows the uni- 
versal scaling S/N ~ E^} in the regime where the three 
structures all have "small-world" property, I ~ \nN. 
Moreover, we have tested other network models [2(| shar- 
ing the property I ~ In TV. The results indicate that the 
"small-world" property plays the key role in the stable 
regime and is responsible for the universal scaling, re- 
gardless of the model details, which can be explained in 
the framework of a mean- field approach [201 ] . 

In summary, we have studied the vibrational thermo- 
dynamic stability of small-world structures. The average 
mean-square displacement S of the structure has been ex- 
pressed as the mean of inverse eigenvalues of its Laplacian 
matrix L. Therefore, the dynamic vibration property is 
closely related to the static structure information. It is 
found that the cross-links suppress S effectively and on 
the small- world network model, there exist two phases: 
an unstable phase where p <C 1/N, S ~ N and a sta- 
ble phase where p > 1/N, S ~ p~ l , i.e., S/N ~ E~*. 
Further, we have tested various data from the PDB and 
find that native proteins belong to the stable phase and 
share the same scaling behavior S/N ~ E~ y . It is be- 
lieved that nature selects proteins under the constraint of 
thermodynamic stability so that proteins can keep their 
specific native fold structure stable for proper function. 
Finally, we have studied S in ER random network model 
and have validated our conjecture that it is the "small- 
world" property that plays a key role in the thermody- 
namic stability of structures and is responsible for the 
universal scaling, S/N ~ E^}, in the stable regime. It 
is also interesting to examine more complex structure 



effects on the thermodynamic stability problem, such as 
scale- free networks [20], hierarchical structures, networks 
with community structure etc. More realistic consider- 
ations such as the effect of random coupling constants, 
anharmonic potentials, or even quantum version of vi- 
bration dynamics are worth further studying. 

The work is supported by the NUS Faculty Research 
Grant No. R-144-000-165-112/101. 



APPENDIX A: DERIVATION OF THE 
CORRELATION MATRIX 

To make this paper self-contained and readable, we 
complement the detailed derivation of Eq. ^ which is 
expressed in terms of G in Fourier transform space. We 
follow Ref. [21] by defining the Fourier transform as 



1 f + °° 

QM = ^ J q(t) ( 

i r +QO 
,M - T J {(*). 



- luJt dt- 



(Al) 
(A2) 



Applying Fourier transform to both sides of Eq. (1), 
one obtains 

- u?MQ = -oLQ - iluTQ + r). (A3) 

Simple algebraic operation yields, 

Q = G- 1 (iLu)r 1 , (A4) 

where matrix G(iu>) = —oj 2 M + iojT + aL as denned in 
text. The two point correlation function is 

(qi(t + T)q k (t)) 



+00 



-oo 



*^{u)rf{J))\G- y {iu)G- x {-iJ)\ lk . 



(A5) 



where * denotes the conjugate transpose and Q*(u>') 
ry*(w')G _1 (— ioj'). Moreover, since 



1 f + °° 
(r)(uj)ri*(u})) = — I <u< 
Zir 



dte~ luJt 



1 



-Hoc 



X _L / dt'e^'mCit')) 



i(u>'-w)t 



k B TT 1 r+°° 
71" 2tt J_ 00 

= **EE_s(w'-w), (A6) 

7T 

substitute Eq. (|A6[) into Eq. (| A5|) and we have: 
(qi{t + T)q k (t)) 

dwe <WT [G -1 (« t ;)rG- 1 Ha;)]ifc. (A7) 



k B T 
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The expression of correlation matrix Eq. ([2]) corresponds 
to the special case r = 0. 



APPENDIX B: INFORMATION OF PROTEINS 
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0.06937 


1CB8 
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28.0511 


1.164 


0.03575 


1A1S 


313 


1088 


21.3477 


1.068 


0.06386 


1HMU 


674 


2341 


21.8915 


0.907 


0.03581 


IADS 


315 


993 


10.7205 


0.5 


0.06807 


1A47 


683 


2350 


13.5361 


0.646 


0.03068 


1A40 


321 


1153 


9.72025 


0.524 


0.05779 


1CDG 


686 


2375 


23.1041 


1.074 


0.03136 


1A54 


321 


1144 


11.6098 


0.601 


0.06018 


1DMT 
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2313 


26.9702 


1.204 
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1A0I 
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1094 


27.2887 


1.109 


0.07412 


1A4G 


780 


2904 


11.4584 


0.591 


0.02486 


3PTE 


347 


1210 


8.13032 


0.366 


0.06402 


1HTY 
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3276 


14.0281 


0.646 


0.02198 
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351 


1117 


34.2736 


1.369 


0.07133 


1KCW 


1017 
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44.0269 


2.13 


0.02032 


1BVW 


360 


1209 


13.0188 


0.652 


0.05547 
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8860 


26.7142 


1.263 


0.00859 


8JDW 


360 


1191 


23.8334 


1.293 


0.05120 


1B0P 


2462 


8936 


6.08348 


0.319 


0.00775 
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1266 


19.8278 


0.859 


0.05964 


1K83 
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55.9118 


2.03 


0.00788 
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0.00813 


1A39 


410 


1474 


21.4742 


1.113 


0.04706 


1150 


3558 


11799 


63.2545 


2.236 


0.00795 



TABLE I: Information of Proteins used in the present study. Size N is 
the number of residues. E cr is the number of cross-links, counted within 
cutoff 7 A. B is the average B-factor over all G a atoms for each protein. 
The estimated ksT /a are collected in Ref. 0. B'/N is the normalized 
B-factor over the size of protein. 
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